Setting Initial Variables

dv = "irt.em.norm"
dv.pos = "pos.irt.em.norm"
dv.pre = "pre.irt.em.norm"
dv.dif = "dif.irt.em.norm"

fatores2 <- c("Gender","Town","Degree","qtl.irt.em.norm")
lfatores2 <- as.list(fatores2)
names(lfatores2) <- fatores2

fatores1 <- c("Group", fatores2)
lfatores1 <- as.list(fatores1)
names(lfatores1) <- fatores1

lfatores <- c(lfatores1)

color <- list()
color[["prepost"]] = c("#ffee65","#f28e2B")
color[["Group"]] = c("#bcbd22","#fd7f6f","#008000")
color[["Gender"]] = c("#FF007F","#4D4DFF")
color[["Town"]] = c("#AA00CC","#00AA99")
color[["Degree"]] = c("#d21820","#f36dff","#aa8882")


level <- list()
level[["Group"]] = c("Ctr.A","Ctr.B","Exp")
level[["Gender"]] = c("Female","Male")
level[["Town"]] = c("Sorocaba - SP","Limoeiro - PE")
level[["Degree"]] = c("3a","4a","5a")
level[["qtl.irt.em.norm"]] = c("1st","2nd","3rd")
level[["qtl.score"]] = c("1st","2nd","3rd")

# ..

color[["Group:Gender"]] = c(
  "Ctr.A:Female"="#ffccbb", "Ctr.B:Female"="#ff99cc", "Exp:Female"="#FF007F",
  "Ctr.A:Male"="#aabbff", "Ctr.B:Male"="#bbaaff", "Exp:Male"="#4D4DFF",
  "Ctr.A.Female"="#ffccbb", "Ctr.B.Female"="#ff99cc", "Exp.Female"="#FF007F",
  "Ctr.A.Male"="#aabbff", "Ctr.B.Male"="#bbaaff", "Exp.Male"="#4D4DFF"
)
color[["Group:Town"]] = c(
  "Ctr.A:Sorocaba - SP"="#AA00FF", "Ctr.B:Sorocaba - SP"="#FF00FF", "Exp:Sorocaba - SP"="#AA00FF",
  "Ctr.A:Limoeiro - PE"="#00EEFF", "Ctr.B:Limoeiro - PE"="#00EECC", "Exp:Limoeiro - PE"="#00CFCF",
  "Ctr.A.Sorocaba - SP"="#AA00FF", "Ctr.B.Sorocaba - SP"="#FF00FF", "Exp.Sorocaba - SP"="#AA00FF",
  "Ctr.A.Limoeiro - PE"="#00EEFF", "Ctr.B.Limoeiro - PE"="#00EECC", "Exp.Limoeiro - PE"="#00CFCF"
)


for (coln in c("vocab")) {
  color[[paste0(coln,".quintile")]] = c("#BF0040","#FF0000","#800080","#0000FF","#4000BF")
  level[[paste0(coln,".quintile")]] = c("1st quintile","2nd quintile","3rd quintile","4th quintile","5th quintile")
  color[[paste0("grupo:",coln,".quintile")]] = c(
    "Experimental.1st quintile"="#BF0040", "Controle.1st quintile"="#d8668c",
    "Experimental.2nd quintile"="#FF0000", "Controle.2nd quintile"="#ff7f7f",
    "Experimental.3rd quintile"="#8fce00", "Controle.3rd quintile"="#ddf0b2",
    "Experimental.4th quintile"="#0000FF", "Controle.4th quintile"="#b2b2ff",
    "Experimental.5th quintile"="#4000BF", "Controle.5th quintile"="#b299e5",
    
    "Experimental:1st quintile"="#BF0040", "Controle:1st quintile"="#d8668c",
    "Experimental:2nd quintile"="#FF0000", "Controle:2nd quintile"="#ff7f7f",
    "Experimental:3rd quintile"="#8fce00", "Controle:3rd quintile"="#ddf0b2",
    "Experimental:4th quintile"="#0000FF", "Controle:4th quintile"="#b2b2ff",
    "Experimental:5th quintile"="#4000BF", "Controle:5th quintile"="#b299e5")
}


gdat <- read_excel("../data/data.xlsx", sheet = "data")
gdat <- gdat[!is.na(gdat[[dv.pre]]) & !is.na(gdat[[dv.pos]]),]

gdat <- gdat[which(gdat$dif.irt.em.norm != 0 & gdat$pre.irt.em.norm != 100),]

dat <- gdat
dat$Group <- factor(dat[["Group"]], level[["Group"]])
for (coln in c(names(lfatores))) {
  if (length(level[[coln]]) > 0)
    plevel = level[[coln]][level[[coln]] %in% unique(dat[[coln]])]
  else
    plevel = unique(dat[[coln]])[!is.na(unique(dat[[coln]]))]
  
  dat[[coln]] <- factor(dat[[coln]], plevel)
}

dat <- dat[,c("ID", names(lfatores), dv.pre, dv.pos, dv.dif)]

dat.long <- rbind(dat, dat)
dat.long$time <- c(rep("pre", nrow(dat)), rep("pos", nrow(dat)))
dat.long$time <- factor(dat.long$time, c("pre","pos"))
dat.long[[dv]] <- c(dat[[dv.pre]], dat[[dv.pos]])


for (f in c("Group", names(lfatores))) {
  if (is.null(color[[f]]) && length(unique(dat[[f]])) > 0) 
      color[[f]] <- distinctColorPalette(length(unique(dat[[f]])))
}

for (f in c(fatores2)) {
  if (is.null(color[[paste0("Group:",f)]]) && length(unique(dat[[f]])) > 0)
    color[[paste0("Group:",f)]] <- distinctColorPalette(
      length(unique(dat[["Group"]]))*length(unique(dat[[f]])))
}

ldat <- list()
laov <- list()
lpwc <- list()
lemms <- list()

Descriptive Statistics of Initial Data

df <- get.descriptives(dat, c(dv.pre, dv.pos, dv.dif), c("Group"),
                       symmetry.test = T, normality.test = F)
df <- plyr::rbind.fill(
  df, do.call(plyr::rbind.fill, lapply(lfatores2, FUN = function(f) {
    if (nrow(dat) > 0 && sum(!is.na(unique(dat[[f]]))) > 1)
      get.descriptives(dat, c(dv.pre,dv.pos), c("Group", f), include.global = F,
                       symmetry.test = T, normality.test = F)
    }))
)
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning in `stats::qt()`:
## ! NaNs produced
## There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning in `stats::qt()`:
## ! NaNs produced
## There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning in `stats::qt()`:
## ! NaNs produced
## There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning in `stats::qt()`:
## ! NaNs produced
df <- df[,c("variable",fatores1[fatores1 %in% colnames(df)],
            colnames(df)[!colnames(df) %in% c(fatores1,"variable")])]
variable Group Gender Town Degree qtl.irt.em.norm n mean median min max sd se ci iqr symmetry skewness kurtosis
pre.irt.em.norm Ctr.A 14 61.082 75.243 23.198 75.243 20.000 5.345 11.547 30.446 NO -0.7087597 -1.3689578
pre.irt.em.norm Ctr.B 68 51.016 47.360 11.805 93.027 22.941 2.782 5.553 41.467 YES -0.1454210 -1.3896233
pre.irt.em.norm Exp 42 50.220 51.925 11.805 75.781 23.735 3.662 7.396 41.467 YES -0.3646594 -1.4015677
pos.irt.em.norm Ctr.A 14 55.233 55.682 0.000 100.000 31.798 8.498 18.359 41.720 YES 0.0127877 -1.2696808
pos.irt.em.norm Ctr.B 68 66.967 68.318 0.000 100.000 32.760 3.973 7.930 58.460 YES -0.4861066 -1.1182797
pos.irt.em.norm Exp 42 66.882 69.560 0.000 100.000 33.243 5.129 10.359 56.467 YES -0.4824511 -1.1061834
dif.irt.em.norm Ctr.A 14 -5.850 -7.947 -48.335 40.593 27.762 7.420 16.029 47.352 YES 0.0723598 -1.3581213
dif.irt.em.norm Ctr.B 68 15.951 24.757 -62.411 88.195 30.115 3.652 7.289 36.834 YES -0.0672311 -0.2059701
dif.irt.em.norm Exp 42 16.662 22.688 -31.710 81.972 25.131 3.878 7.831 33.925 YES 0.2342158 -0.3658467
pre.irt.em.norm Ctr.A Female 9 61.470 75.243 23.198 75.243 21.354 7.118 16.414 30.446 NO -0.7560725 -1.3824688
pre.irt.em.norm Ctr.A Male 5 60.385 72.864 33.776 75.243 19.674 8.798 24.428 30.446 YES -0.3709082 -2.0895003
pre.irt.em.norm Ctr.B Female 36 49.292 44.797 11.805 75.243 23.402 3.900 7.918 41.467 YES -0.1219004 -1.5262112
pre.irt.em.norm Ctr.B Male 31 53.915 49.610 11.805 93.027 22.326 4.010 8.189 41.467 YES -0.2119987 -1.2883002
pre.irt.em.norm Ctr.B 1 23.198 23.198 23.198 23.198 0.000 few data 0.0000000 0.0000000
pre.irt.em.norm Exp Female 18 51.025 46.078 11.805 75.243 23.308 5.494 11.591 41.467 YES -0.2208823 -1.6169865
pre.irt.em.norm Exp Male 20 51.809 60.672 11.805 75.781 24.233 5.419 11.342 38.651 NO -0.6328965 -1.2099954
pre.irt.em.norm Exp 4 38.650 33.776 11.805 75.243 26.503 13.251 42.172 15.860 few data 0.0000000 0.0000000
pos.irt.em.norm Ctr.A Female 9 59.915 63.791 0.000 100.000 32.883 10.961 25.276 23.494 YES -0.4341981 -1.0875306
pos.irt.em.norm Ctr.A Male 5 46.805 43.533 22.432 100.000 31.390 14.038 38.976 19.005 NO 0.8077339 -1.2271576
pos.irt.em.norm Ctr.B Female 36 60.168 68.318 0.000 100.000 35.146 5.858 11.892 77.790 YES -0.1804006 -1.5343996
pos.irt.em.norm Ctr.B Male 31 77.024 70.801 22.210 100.000 25.421 4.566 9.324 36.209 NO -0.5950995 -0.9043817
pos.irt.em.norm Ctr.B 1 0.000 0.000 0.000 0.000 0.000 few data 0.0000000 0.0000000
pos.irt.em.norm Exp Female 18 68.564 68.318 22.210 100.000 27.260 6.425 13.556 55.510 YES -0.1592713 -1.4519300
pos.irt.em.norm Exp Male 20 71.346 87.741 0.000 100.000 34.888 7.801 16.328 50.953 NO -0.7836514 -0.7802592
pos.irt.em.norm Exp 4 36.989 22.210 3.538 100.000 42.919 21.460 68.294 24.115 few data 0.0000000 0.0000000
pre.irt.em.norm Ctr.A Limoeiro - PE 14 61.082 75.243 23.198 75.243 20.000 5.345 11.547 30.446 NO -0.7087597 -1.3689578
pre.irt.em.norm Ctr.B Sorocaba - SP 26 59.953 70.801 28.443 75.243 17.343 3.401 7.005 27.883 YES -0.4899631 -1.5257343
pre.irt.em.norm Ctr.B Limoeiro - PE 42 45.484 35.424 11.805 93.027 24.385 3.763 7.599 49.400 YES 0.2276511 -1.3824847
pre.irt.em.norm Exp Sorocaba - SP 13 61.318 68.318 44.797 75.243 13.163 3.651 7.954 23.706 YES -0.2828356 -1.8876353
pre.irt.em.norm Exp Limoeiro - PE 29 45.245 35.152 11.805 75.781 25.843 4.799 9.830 52.453 YES -0.0014812 -1.7087455
pos.irt.em.norm Ctr.A Limoeiro - PE 14 55.233 55.682 0.000 100.000 31.798 8.498 18.359 41.720 YES 0.0127877 -1.2696808
pos.irt.em.norm Ctr.B Sorocaba - SP 26 89.615 100.000 22.210 100.000 19.487 3.822 7.871 18.568 NO -1.8262078 2.9780354
pos.irt.em.norm Ctr.B Limoeiro - PE 42 52.947 55.575 0.000 100.000 31.553 4.869 9.833 48.591 YES 0.0566616 -1.1944851
pos.irt.em.norm Exp Sorocaba - SP 13 90.201 100.000 29.058 100.000 21.216 5.884 12.820 0.000 NO -1.8769367 2.4478110
pos.irt.em.norm Exp Limoeiro - PE 29 56.428 49.610 0.000 100.000 32.551 6.045 12.382 41.822 YES -0.0905270 -1.1735271
pre.irt.em.norm Ctr.A 3a 2 75.243 75.243 75.243 75.243 0.000 0.000 0.000 0.000 few data 0.0000000 0.0000000
pre.irt.em.norm Ctr.A 4a 5 50.451 44.797 23.198 75.243 23.885 10.682 29.657 41.467 YES 0.0794730 -2.1486168
pre.irt.em.norm Ctr.A 5a 7 64.630 75.243 33.776 75.243 17.624 6.661 16.299 16.412 NO -0.8307167 -1.3470817
pre.irt.em.norm Ctr.B 3a 28 45.038 33.776 11.805 93.027 27.246 5.149 10.565 52.045 YES 0.1962780 -1.6228339
pre.irt.em.norm Ctr.B 4a 21 54.257 44.797 23.198 75.243 19.668 4.292 8.953 41.467 YES -0.0272718 -1.8255975
pre.irt.em.norm Ctr.B 5a 19 56.244 49.610 23.198 75.243 17.860 4.097 8.608 27.975 YES -0.2720853 -1.4775977
pre.irt.em.norm Exp 3a 19 46.374 44.797 11.805 75.781 25.261 5.795 12.176 46.960 YES -0.1079885 -1.6518026
pre.irt.em.norm Exp 4a 15 55.156 63.791 18.028 75.243 19.498 5.034 10.798 26.136 NO -0.5340963 -1.2192908
pre.irt.em.norm Exp 5a 8 50.097 58.831 11.805 75.243 28.391 10.038 23.735 46.960 YES -0.3129260 -1.8764432
pos.irt.em.norm Ctr.A 3a 2 100.000 100.000 100.000 100.000 0.000 0.000 0.000 0.000 few data 0.0000000 0.0000000
pos.irt.em.norm Ctr.A 4a 5 50.393 43.533 22.210 100.000 32.648 14.601 40.538 41.360 YES 0.4592566 -1.7092278
pos.irt.em.norm Ctr.A 5a 7 45.899 47.572 0.000 71.066 26.285 9.935 24.310 33.265 NO -0.5455584 -1.3435247
pos.irt.em.norm Ctr.B 3a 28 46.955 43.533 0.000 100.000 31.513 5.955 12.219 48.591 YES 0.2482926 -1.1163764
pos.irt.em.norm Ctr.B 4a 21 82.305 100.000 22.210 100.000 23.282 5.081 10.598 31.682 NO -0.9366991 -0.1662967
pos.irt.em.norm Ctr.B 5a 19 79.506 100.000 22.210 100.000 28.972 6.647 13.964 31.682 NO -1.0339818 -0.4318475
pos.irt.em.norm Exp 3a 19 54.597 49.610 0.000 100.000 32.727 7.508 15.774 45.207 YES -0.0652619 -1.2265291
pos.irt.em.norm Exp 4a 15 80.445 100.000 22.210 100.000 28.230 7.289 15.633 31.682 NO -0.9479436 -0.7231581
pos.irt.em.norm Exp 5a 8 70.625 84.159 3.538 100.000 36.121 12.771 30.198 51.909 NO -0.6092978 -1.2400616
pre.irt.em.norm Ctr.A 1st 1 23.198 23.198 23.198 23.198 0.000 few data 0.0000000 0.0000000
pre.irt.em.norm Ctr.A 2nd 5 46.002 44.797 33.776 72.864 15.995 7.153 19.861 11.021 few data 0.0000000 0.0000000
pre.irt.em.norm Ctr.A 3rd 8 75.243 75.243 75.243 75.243 0.000 0.000 0.000 0.000 few data 0.0000000 0.0000000
pre.irt.em.norm Ctr.B 1st 12 17.939 17.501 11.805 28.443 6.567 1.896 4.173 11.393 few data 0.0000000 0.0000000
pre.irt.em.norm Ctr.B 2nd 33 45.620 44.797 33.776 72.864 13.715 2.387 4.863 13.583 NO 0.9057667 -0.6149984
pre.irt.em.norm Ctr.B 3rd 23 76.016 75.243 75.243 93.027 3.708 0.773 1.604 0.000 few data 0.0000000 0.0000000
pre.irt.em.norm Exp 1st 9 14.983 11.805 11.805 23.198 4.979 1.660 3.827 6.223 NO 0.7796101 -1.3784196
pre.irt.em.norm Exp 2nd 21 50.997 47.360 33.776 72.864 14.441 3.151 6.573 28.639 YES 0.1626342 -1.5980329
pre.irt.em.norm Exp 3rd 12 75.288 75.243 75.243 75.781 0.155 0.045 0.099 0.000 few data 0.0000000 0.0000000
pos.irt.em.norm Ctr.A 1st 1 63.791 63.791 63.791 63.791 0.000 few data 0.0000000 0.0000000
pos.irt.em.norm Ctr.A 2nd 5 22.541 22.432 0.000 43.533 15.434 6.902 19.164 2.319 YES -0.1193175 -1.4055909
pos.irt.em.norm Ctr.A 3rd 8 74.595 70.933 43.533 100.000 23.222 8.210 19.414 40.264 YES -0.0187306 -1.8274973
pos.irt.em.norm Ctr.B 1st 12 37.670 22.210 0.000 100.000 35.017 10.109 22.249 29.251 NO 0.7927758 -0.9741981
pos.irt.em.norm Ctr.B 2nd 33 69.806 68.318 0.000 100.000 30.844 5.369 10.937 56.467 YES -0.4999557 -1.0393866
pos.irt.em.norm Ctr.B 3rd 23 78.180 70.801 12.832 100.000 25.685 5.356 11.107 31.682 NO -0.9746332 0.1539423
pos.irt.em.norm Exp 1st 9 35.637 22.210 0.000 100.000 36.024 12.008 27.691 64.780 NO 0.5147071 -1.4045460
pos.irt.em.norm Exp 2nd 21 70.943 70.801 22.210 100.000 29.244 6.382 13.312 52.640 YES -0.3005251 -1.5268078
pos.irt.em.norm Exp 3rd 12 83.207 100.000 43.533 100.000 22.343 6.450 14.196 32.813 NO -0.6458141 -1.3681401

One-way factor analysis for: irt.em.norm ~ Group

pdat = remove_group_data(dat[!is.na(dat[["Group"]]),], "dif.irt.em.norm", "Group")

pdat.long <- rbind(pdat[,c("ID","Group")], pdat[,c("ID","Group")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["irt.em.norm"]] <- c(pdat[["pre.irt.em.norm"]], pdat[["pos.irt.em.norm"]])

y.position.min <- abs(
  max(pdat.long[["irt.em.norm"]])
  - min(pdat.long[["irt.em.norm"]]))/20

lvars = as.list(c("dif.irt.em.norm","pos.irt.em.norm","pre.irt.em.norm"))
names(lvars) = unlist(lvars)

Pre-test and Post-test PairWise comparisons for: irt.em.norm ~ Group

pwc.long <- group_by(pdat.long, Group) %>%
  pairwise_wilcox_test(irt.em.norm ~ time, detailed = T)

df <- pwc.long[,c(".y.","Group","group1","group2","n1","n2","estimate",
                  "statistic","p.adj","p.adj.signif")]
.y. Group group1 group2 n1 n2 estimate statistic p.adj p.adj.signif
irt.em.norm Ctr.A pre pos 14 14 7.672307 122.0 0.274 ns
irt.em.norm Ctr.B pre pos 68 68 -24.757101 1734.5 0.011 *
irt.em.norm Exp pre pos 42 42 -24.757119 617.0 0.017 * 
stat.test <- pwc.long %>% add_xy_position(x = "Group", fun = "max")
stat.test$y.position <- stat.test$y.position + y.position.min

ggboxplot(pdat.long, x = "Group", y = "irt.em.norm",
          palette = color$prepost, fill = "time") +
  stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = T,
                     label = "{ p.adj } ({ p.adj.signif })") + xlab("")

stat.test <- pwc.long %>% add_xy_position(x = "time", fun = "max")
stat.test$y.position <- stat.test$y.position + y.position.min

gg <- ggline(
  pdat.long, x = "time", y = "irt.em.norm", size = 1.5,
  facet.by = "Group", add = c("mean_se"), color = "Group",
  position = position_dodge(width = 0.3), palette = color[["Group"]])

pdat.long$xj = jitter(as.numeric(pdat.long[["time"]]), amount=.1)
pdat.long$yj = jitter(pdat.long[["irt.em.norm"]], amount = .01)

gg + geom_point(
  data = pdat.long, aes_string(x="xj",y="yj", color = "Group"), size=0.5) +
  stat_pvalue_manual(
    stat.test, tip.length = 0, hide.ns = T, label.size = 5,
    position = position_dodge(width = 0.3),
    label = "{ p.adj } ({ p.adj.signif })") + xlab("") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

Kruskal and Wilcoxon PairWise comparisons for: irt.em.norm ~ Group

kt <- lapply(lvars, FUN = function(x) {
  kruskal_test(pdat, as.formula(paste0(x," ~ Group")))  
})

df <- do.call(rbind.fill, lapply(lvars, function(x) {
  add_significance(merge(
    kt[[x]], kruskal_effsize(pdat, as.formula(paste0(x," ~ Group"))),
    by = c(".y.","n"), suffixes = c("",".ez")))
}))

df <- df[,c(".y.","n","df","statistic","p","p.signif","effsize","magnitude")]
.y. n df statistic p p.signif effsize magnitude
dif.irt.em.norm 124 2 6.413810 0.0405 * 0.0364778 small
pos.irt.em.norm 124 2 1.407533 0.4950 ns -0.0048964 small
pre.irt.em.norm 124 2 2.832761 0.2430 ns 0.0068823 small 
pwc <- lapply(lvars, FUN = function(x) {
  pairwise_wilcox_test(pdat, as.formula(paste0(x," ~ Group")))  
})

df <- do.call(rbind.fill, pwc)
.y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
dif.irt.em.norm Ctr.A Ctr.B 14 68 280.0 0.016 0.047 *
dif.irt.em.norm Ctr.A Exp 14 42 170.5 0.020 0.047 *
dif.irt.em.norm Ctr.B Exp 68 42 1453.0 0.880 0.880 ns
pos.irt.em.norm Ctr.A Ctr.B 14 68 392.0 0.290 0.735 ns
pos.irt.em.norm Ctr.A Exp 14 42 233.5 0.245 0.735 ns
pos.irt.em.norm Ctr.B Exp 68 42 1405.0 0.886 0.886 ns
pre.irt.em.norm Ctr.A Ctr.B 14 68 597.5 0.125 0.306 ns
pre.irt.em.norm Ctr.A Exp 14 42 379.0 0.102 0.306 ns
pre.irt.em.norm Ctr.B Exp 68 42 1464.5 0.822 0.822 ns 
plots <- lapply(lvars, FUN = function(y) {
  stat.test <- pwc[[y]] %>% add_xy_position(x = "Group")
  stat.test$y.position <- stat.test$y.position + y.position.min
  ggboxplot(pdat, x = "Group", y = y, fill = "Group",
            palette = color[["Group"]]) +
    stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = T, label.size = 5,
                       label="{ p.adj } ({ p.adj.signif })") + xlab("")
})
egg::ggarrange(plots[["pre.irt.em.norm"]], plots[["pos.irt.em.norm"]], nrow = 1)

plots[["dif.irt.em.norm"]] +
  labs(subtitle = get_test_label(kt[["dif.irt.em.norm"]], detailed = T),
       caption = get_pwc_label(pwc[["dif.irt.em.norm"]])) +
  ylab("irt.em.norm (dif)")  +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

Two-way factor analysis for: irt.em.norm ~ Group:Gender

pdat = remove_group_data(
  dat[!is.na(dat[["Group"]]) & !is.na(dat[["Gender"]]),],
  "dif.irt.em.norm", c("Group","Gender"))

pdat.long <- rbind(pdat[,c("ID","Group","Gender")],
                   pdat[,c("ID","Group","Gender")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["irt.em.norm"]] <- c(pdat[["pre.irt.em.norm"]], pdat[["pos.irt.em.norm"]])

y.position.min <- abs(
  max(pdat.long[["irt.em.norm"]])
  - min(pdat.long[["irt.em.norm"]]))/20

lvars = as.list(c("dif.irt.em.norm","pos.irt.em.norm","pre.irt.em.norm"))
names(lvars) = unlist(lvars)

Pre-test and Post-test PairWise comparisons for: irt.em.norm ~ Group:Gender

pwc.long <- group_by(pdat.long, Group:Gender) %>%
  pairwise_wilcox_test(irt.em.norm ~ time, detailed = T)

df <- pwc.long[,c(".y.","Group:Gender","group1","group2","n1","n2","estimate",
                  "statistic","p.adj","p.adj.signif")]
.y. Group:Gender group1 group2 n1 n2 estimate statistic p.adj p.adj.signif
irt.em.norm Ctr.A:Female pre pos 9 9 4.441848 48.0 0.529 ns
irt.em.norm Ctr.A:Male pre pos 5 5 22.365539 18.0 0.293 ns
irt.em.norm Ctr.B:Female pre pos 36 36 -10.405021 574.5 0.408 ns
irt.em.norm Ctr.B:Male pre pos 31 31 -24.757179 270.0 0.003 **
irt.em.norm Exp:Female pre pos 18 18 -23.521266 119.0 0.176 ns
irt.em.norm Exp:Male pre pos 20 20 -24.757198 120.5 0.031 * 
stat.test <- pwc.long %>% add_xy_position(x = "Group:Gender", fun = "max")
sidx = which(stat.test$p.adj.signif != "ns")
stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))

pdat.long[[paste0(c("Group","Gender"), collapse = ":")]] = apply(
  pdat.long[, c("Group","Gender")], 1, paste0, collapse = ":")

ggboxplot(pdat.long, x = "Group:Gender", y = "irt.em.norm",
          palette = color$prepost, fill = "time") +
  stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = T,
                     label = "{ p.adj } ({ p.adj.signif })") + xlab("")

pwc.long <- group_by(pdat.long, Group, Gender) %>%
  pairwise_wilcox_test(irt.em.norm ~ time, detailed = T)

stat.test <- pwc.long %>% add_xy_position(x = "time", fun = "mean_se")
sidx = which(stat.test$p.adj.signif != "ns")
stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))

gg <- ggline(
  pdat.long, x = "time", y = "irt.em.norm",
  color = "Gender", linetype = "Gender", shape = "Gender", size = 1.5,
  facet.by = "Group", add = c("mean_se"),
  position = position_dodge(width = 0.3), palette = color[["Gender"]])

pdat.long$xj = jitter(as.numeric(pdat.long[["time"]]), amount=.1)
pdat.long$yj = jitter(pdat.long[["irt.em.norm"]], amount = .01)

gg + geom_point(
  data = pdat.long, aes_string(x="xj",y="yj",colour="Gender"), size=0.5) +
  stat_pvalue_manual(
    stat.test, tip.length = 0, hide.ns = T, label.size = 5,
    position = position_dodge(width = 0.3), color = "Gender",
    label = "{ p.adj } ({ p.adj.signif })") + xlab("") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

Scheirer and Wilcoxon PairWise comparisons for: irt.em.norm ~ Group:Gender

sch <- lapply(lvars, FUN = function(x) {
  scheirer.test(pdat, x, c("Group","Gender"), as.table = T) 
})
df <- do.call(rbind.fill, sch)
var Effect Df Sum Sq H p.value p.value.signif
dif.irt.em.norm Group 2 7658.5419 6.4692197 0.0393756 *
dif.irt.em.norm Gender 1 1804.5449 1.5243107 0.2169683 ns
dif.irt.em.norm Group:Gender 2 1691.0919 1.4284762 0.4895650 ns
dif.irt.em.norm Residuals 113 127829.9913
pos.irt.em.norm Group 2 1916.2193 1.7152319 0.4241721 ns
pos.irt.em.norm Gender 1 2514.9765 2.2511870 0.1335120 ns
pos.irt.em.norm Group:Gender 2 2270.3878 2.0322526 0.3619945 ns
pos.irt.em.norm Residuals 113 124655.0630
pre.irt.em.norm Group 2 2942.0016 2.5847863 0.2746128 ns
pre.irt.em.norm Gender 1 289.4244 0.2542827 0.6140755 ns
pre.irt.em.norm Group:Gender 2 465.9460 0.4093712 0.8149035 ns
pre.irt.em.norm Residuals 113 130745.8694 
pwc <- lapply(lvars, FUN = function(x) {
  list(
    Group = tryCatch(pairwise_wilcox_test(group_by(pdat, Gender),
                                 as.formula(paste0(x," ~ Group")))
                         , error = function(e) NULL),
    Gender = tryCatch(pairwise_wilcox_test(group_by(pdat, Group),
                                 as.formula(paste0(x," ~ Gender")))
                         , error = function(e) NULL)
  )
})

df <- do.call(rbind.fill, lapply(pwc, FUN =  function(x) {
  do.call(rbind.fill, x)
}))

ivs = c()
if ("Group" %in% colnames(df)) ivs = c(ivs, "Group")
if ("Gender" %in% colnames(df)) ivs = c(ivs, "Gender")
df <- df[,c(".y.",ivs,"group1","group2","n1","n2",
            "statistic","p.adj","p.adj.signif")]
.y. Group Gender group1 group2 n1 n2 statistic p.adj p.adj.signif
dif.irt.em.norm Female Ctr.A Ctr.B 9 36 118.5 0.489 ns
dif.irt.em.norm Female Ctr.A Exp 9 18 53.5 0.489 ns
dif.irt.em.norm Female Ctr.B Exp 36 18 308.0 0.775 ns
dif.irt.em.norm Male Ctr.A Ctr.B 5 31 26.5 0.062 ns
dif.irt.em.norm Male Ctr.A Exp 5 20 18.0 0.064 ns
dif.irt.em.norm Male Ctr.B Exp 31 20 323.5 0.802 ns
dif.irt.em.norm Ctr.A Female Male 9 5 29.0 0.422 ns
dif.irt.em.norm Ctr.B Female Male 36 31 446.5 0.162 ns
dif.irt.em.norm Exp Female Male 18 20 157.0 0.510 ns
pos.irt.em.norm Female Ctr.A Ctr.B 9 36 159.0 1.000 ns
pos.irt.em.norm Female Ctr.A Exp 9 18 71.0 1.000 ns
pos.irt.em.norm Female Ctr.B Exp 36 18 279.0 1.000 ns
pos.irt.em.norm Male Ctr.A Ctr.B 5 31 35.5 0.139 ns
pos.irt.em.norm Male Ctr.A Exp 5 20 28.0 0.252 ns
pos.irt.em.norm Male Ctr.B Exp 31 20 319.5 0.853 ns
pos.irt.em.norm Ctr.A Female Male 9 5 29.0 0.421 ns
pos.irt.em.norm Ctr.B Female Male 36 31 421.5 0.076 ns
pos.irt.em.norm Exp Female Male 18 20 157.0 0.494 ns
pre.irt.em.norm Female Ctr.A Ctr.B 9 36 213.5 0.399 ns
pre.irt.em.norm Female Ctr.A Exp 9 18 106.0 0.399 ns
pre.irt.em.norm Female Ctr.B Exp 36 18 313.0 0.844 ns
pre.irt.em.norm Male Ctr.A Ctr.B 5 31 88.0 1.000 ns
pre.irt.em.norm Male Ctr.A Exp 5 20 59.5 1.000 ns
pre.irt.em.norm Male Ctr.B Exp 31 20 323.0 1.000 ns
pre.irt.em.norm Ctr.A Female Male 9 5 26.0 0.657 ns
pre.irt.em.norm Ctr.B Female Male 36 31 497.5 0.440 ns
pre.irt.em.norm Exp Female Male 18 20 180.0 1.000 ns 
plots <- lapply(lvars, FUN = function(y) {
  livs = list("Group", "Gender")
  names(livs) = unlist(livs)
  lapply(livs, FUN = function(x) {
    iv2 = setdiff(names(livs), x)
    if (!is.null(pwc[[y]][[iv2]])) {
      stat.test <- pwc[[y]][[iv2]] %>% add_xy_position(x = x, fun = "max")
      sidx = which(stat.test$p.adj.signif != "ns")
      stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))
      
      ggboxplot(pdat, x = x, y = y, fill = iv2, palette = color[[iv2]]) +
        stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = T, label.size = 5,
                           label="{ p.adj } ({ p.adj.signif })") + xlab("")
    }
  })
})
if (!is.null(plots[["pre.irt.em.norm"]][["Group"]]) &&
    !is.null(plots[["pos.irt.em.norm"]][["Group"]])) {
  egg::ggarrange(plots[["pre.irt.em.norm"]][["Group"]],
                 plots[["pos.irt.em.norm"]][["Group"]], nrow = 1)  
}

if (!is.null(plots[["pre.irt.em.norm"]][["Gender"]]) &&
    !is.null(plots[["pos.irt.em.norm"]][["Gender"]])) {
  egg::ggarrange(plots[["pre.irt.em.norm"]][["Gender"]],
                 plots[["pos.irt.em.norm"]][["Gender"]], nrow = 1)
}

psch = sch[["dif.irt.em.norm"]]
idx = which(psch$Effect == "Group:Gender") 

dof = floor(as.double(psch$Df[idx]))
dof.res = floor(as.double(psch$Df[which(psch$Effect == "Residuals")]))
statistic = round(as.double(psch$H[idx]), 3)
p = round(as.double(psch[["p.value"]][idx]), 3)
pval = ifelse(p < 0.001,paste0(" , p<0.001"),paste0(" , p=",p))

if (!is.null(plots[["dif.irt.em.norm"]][["Group"]]))
  plots[["dif.irt.em.norm"]][["Group"]] +
    labs(subtitle = paste0("Scheirer-Ray-Hare H(", dof, ",", 
          dof.res, ")=", statistic, pval),
         caption = get_pwc_label(pwc[["dif.irt.em.norm"]][["Gender"]])) +
    ylab("irt.em.norm (dif)") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

psch = sch[["dif.irt.em.norm"]]
idx = which(psch$Effect == "Group:Gender") 

dof = floor(as.double(psch$Df[idx]))
dof.res = floor(as.double(psch$Df[which(psch$Effect == "Residuals")]))
statistic = round(as.double(psch$H[idx]), 3)
p = round(as.double(psch[["p.value"]][idx]), 3)
pval = ifelse(p < 0.001,paste0(" , p<0.001"),paste0(" , p=",p))

if (!is.null(plots[["dif.irt.em.norm"]][["Gender"]]))
  plots[["dif.irt.em.norm"]][["Gender"]] +
    labs(subtitle = paste0("Scheirer-Ray-Hare H(", dof, ",", 
          dof.res, ")=", statistic, pval),
         caption = get_pwc_label(pwc[["dif.irt.em.norm"]][["Group"]])) +
    ylab("irt.em.norm (dif)") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

Two-way factor analysis for: irt.em.norm ~ Group:Town

pdat = remove_group_data(
  dat[!is.na(dat[["Group"]]) & !is.na(dat[["Town"]]),],
  "dif.irt.em.norm", c("Group","Town"))

pdat.long <- rbind(pdat[,c("ID","Group","Town")],
                   pdat[,c("ID","Group","Town")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["irt.em.norm"]] <- c(pdat[["pre.irt.em.norm"]], pdat[["pos.irt.em.norm"]])

y.position.min <- abs(
  max(pdat.long[["irt.em.norm"]])
  - min(pdat.long[["irt.em.norm"]]))/20

lvars = as.list(c("dif.irt.em.norm","pos.irt.em.norm","pre.irt.em.norm"))
names(lvars) = unlist(lvars)

Pre-test and Post-test PairWise comparisons for: irt.em.norm ~ Group:Town

pwc.long <- group_by(pdat.long, Group:Town) %>%
  pairwise_wilcox_test(irt.em.norm ~ time, detailed = T)

df <- pwc.long[,c(".y.","Group:Town","group1","group2","n1","n2","estimate",
                  "statistic","p.adj","p.adj.signif")]
.y. Group:Town group1 group2 n1 n2 estimate statistic p.adj p.adj.signif
irt.em.norm Ctr.A:Limoeiro - PE pre pos 14 14 7.672307 122.0 0.274000 ns
irt.em.norm Ctr.B:Sorocaba - SP pre pos 26 26 -26.718532 104.5 0.000011 ****
irt.em.norm Ctr.B:Limoeiro - PE pre pos 42 42 -7.763449 829.0 0.637000 ns
irt.em.norm Exp:Sorocaba - SP pre pos 13 13 -29.199082 20.5 0.000768 ***
irt.em.norm Exp:Limoeiro - PE pre pos 29 29 -10.461074 350.0 0.275000 ns 
stat.test <- pwc.long %>% add_xy_position(x = "Group:Town", fun = "max")
sidx = which(stat.test$p.adj.signif != "ns")
stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))

pdat.long[[paste0(c("Group","Town"), collapse = ":")]] = apply(
  pdat.long[, c("Group","Town")], 1, paste0, collapse = ":")

ggboxplot(pdat.long, x = "Group:Town", y = "irt.em.norm",
          palette = color$prepost, fill = "time") +
  stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = T,
                     label = "{ p.adj } ({ p.adj.signif })") + xlab("")

pwc.long <- group_by(pdat.long, Group, Town) %>%
  pairwise_wilcox_test(irt.em.norm ~ time, detailed = T)

stat.test <- pwc.long %>% add_xy_position(x = "time", fun = "mean_se")
sidx = which(stat.test$p.adj.signif != "ns")
stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))

gg <- ggline(
  pdat.long, x = "time", y = "irt.em.norm",
  color = "Town", linetype = "Town", shape = "Town", size = 1.5,
  facet.by = "Group", add = c("mean_se"),
  position = position_dodge(width = 0.3), palette = color[["Town"]])

pdat.long$xj = jitter(as.numeric(pdat.long[["time"]]), amount=.1)
pdat.long$yj = jitter(pdat.long[["irt.em.norm"]], amount = .01)

gg + geom_point(
  data = pdat.long, aes_string(x="xj",y="yj",colour="Town"), size=0.5) +
  stat_pvalue_manual(
    stat.test, tip.length = 0, hide.ns = T, label.size = 5,
    position = position_dodge(width = 0.3), color = "Town",
    label = "{ p.adj } ({ p.adj.signif })") + xlab("") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

Scheirer and Wilcoxon PairWise comparisons for: irt.em.norm ~ Group:Town

sch <- lapply(lvars, FUN = function(x) {
  scheirer.test(pdat, x, c("Group","Town"), as.table = T) 
})
df <- do.call(rbind.fill, sch)
var Effect Df Sum Sq H p.value p.value.signif
dif.irt.em.norm Group 2 2.958046e+03 2.3022728 0.3162771 ns
dif.irt.em.norm Town 1 1.968526e+04 15.3212058 0.0000907 ****
dif.irt.em.norm Group:Town 1 3.822288e+00 0.0029749 0.9565027 ns
dif.irt.em.norm Residuals 119 1.301052e+05
pos.irt.em.norm Group 2 4.506456e+02 0.3704642 0.8309114 ns
pos.irt.em.norm Town 1 3.811459e+04 31.3330262 0.0000000 ****
pos.irt.em.norm Group:Town 1 4.492732e+01 0.0369336 0.8476003 ns
pos.irt.em.norm Residuals 119 1.097498e+05
pre.irt.em.norm Group 2 6.565912e+03 5.3087169 0.0703440 ns
pre.irt.em.norm Town 1 9.930455e+03 8.0290411 0.0046033 **
pre.irt.em.norm Group:Town 1 1.459214e+01 0.0117981 0.9135045 ns
pre.irt.em.norm Residuals 119 1.386798e+05 
pwc <- lapply(lvars, FUN = function(x) {
  list(
    Group = tryCatch(pairwise_wilcox_test(group_by(pdat, Town),
                                 as.formula(paste0(x," ~ Group")))
                         , error = function(e) NULL),
    Town = tryCatch(pairwise_wilcox_test(group_by(pdat, Group),
                                 as.formula(paste0(x," ~ Town")))
                         , error = function(e) NULL)
  )
})

df <- do.call(rbind.fill, lapply(pwc, FUN =  function(x) {
  do.call(rbind.fill, x)
}))

ivs = c()
if ("Group" %in% colnames(df)) ivs = c(ivs, "Group")
if ("Town" %in% colnames(df)) ivs = c(ivs, "Town")
df <- df[,c(".y.",ivs,"group1","group2","n1","n2",
            "statistic","p.adj","p.adj.signif")]
.y. Town group1 group2 n1 n2 statistic p.adj p.adj.signif
dif.irt.em.norm Sorocaba - SP Ctr.B Exp 26 13 171.5 0.952 ns
dif.irt.em.norm Limoeiro - PE Ctr.A Ctr.B 14 42 220.0 0.330 ns
dif.irt.em.norm Limoeiro - PE Ctr.A Exp 14 29 141.0 0.330 ns
dif.irt.em.norm Limoeiro - PE Ctr.B Exp 42 29 575.0 0.695 ns
pos.irt.em.norm Sorocaba - SP Ctr.B Exp 26 13 161.0 0.771 ns
pos.irt.em.norm Limoeiro - PE Ctr.A Ctr.B 14 42 321.0 1.000 ns
pos.irt.em.norm Limoeiro - PE Ctr.A Exp 14 29 200.5 1.000 ns
pos.irt.em.norm Limoeiro - PE Ctr.B Exp 42 29 559.0 1.000 ns
pre.irt.em.norm Sorocaba - SP Ctr.B Exp 26 13 181.5 0.715 ns
pre.irt.em.norm Limoeiro - PE Ctr.A Ctr.B 14 42 398.5 0.127 ns
pre.irt.em.norm Limoeiro - PE Ctr.A Exp 14 29 274.0 0.127 ns
pre.irt.em.norm Limoeiro - PE Ctr.B Exp 42 29 621.5 0.886 ns 
plots <- lapply(lvars, FUN = function(y) {
  livs = list("Group", "Town")
  names(livs) = unlist(livs)
  lapply(livs, FUN = function(x) {
    iv2 = setdiff(names(livs), x)
    if (!is.null(pwc[[y]][[iv2]])) {
      stat.test <- pwc[[y]][[iv2]] %>% add_xy_position(x = x, fun = "max")
      sidx = which(stat.test$p.adj.signif != "ns")
      stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))
      
      ggboxplot(pdat, x = x, y = y, fill = iv2, palette = color[[iv2]]) +
        stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = T, label.size = 5,
                           label="{ p.adj } ({ p.adj.signif })") + xlab("")
    }
  })
})
if (!is.null(plots[["pre.irt.em.norm"]][["Group"]]) &&
    !is.null(plots[["pos.irt.em.norm"]][["Group"]])) {
  egg::ggarrange(plots[["pre.irt.em.norm"]][["Group"]],
                 plots[["pos.irt.em.norm"]][["Group"]], nrow = 1)  
}
if (!is.null(plots[["pre.irt.em.norm"]][["Town"]]) &&
    !is.null(plots[["pos.irt.em.norm"]][["Town"]])) {
  egg::ggarrange(plots[["pre.irt.em.norm"]][["Town"]],
                 plots[["pos.irt.em.norm"]][["Town"]], nrow = 1)
}

psch = sch[["dif.irt.em.norm"]]
idx = which(psch$Effect == "Group:Town") 

dof = floor(as.double(psch$Df[idx]))
dof.res = floor(as.double(psch$Df[which(psch$Effect == "Residuals")]))
statistic = round(as.double(psch$H[idx]), 3)
p = round(as.double(psch[["p.value"]][idx]), 3)
pval = ifelse(p < 0.001,paste0(" , p<0.001"),paste0(" , p=",p))

if (!is.null(plots[["dif.irt.em.norm"]][["Group"]]))
  plots[["dif.irt.em.norm"]][["Group"]] +
    labs(subtitle = paste0("Scheirer-Ray-Hare H(", dof, ",", 
          dof.res, ")=", statistic, pval),
         caption = get_pwc_label(pwc[["dif.irt.em.norm"]][["Town"]])) +
    ylab("irt.em.norm (dif)") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))
psch = sch[["dif.irt.em.norm"]]
idx = which(psch$Effect == "Group:Town") 

dof = floor(as.double(psch$Df[idx]))
dof.res = floor(as.double(psch$Df[which(psch$Effect == "Residuals")]))
statistic = round(as.double(psch$H[idx]), 3)
p = round(as.double(psch[["p.value"]][idx]), 3)
pval = ifelse(p < 0.001,paste0(" , p<0.001"),paste0(" , p=",p))

if (!is.null(plots[["dif.irt.em.norm"]][["Town"]]))
  plots[["dif.irt.em.norm"]][["Town"]] +
    labs(subtitle = paste0("Scheirer-Ray-Hare H(", dof, ",", 
          dof.res, ")=", statistic, pval),
         caption = get_pwc_label(pwc[["dif.irt.em.norm"]][["Group"]])) +
    ylab("irt.em.norm (dif)") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

Two-way factor analysis for: irt.em.norm ~ Group:Degree

pdat = remove_group_data(
  dat[!is.na(dat[["Group"]]) & !is.na(dat[["Degree"]]),],
  "dif.irt.em.norm", c("Group","Degree"))

pdat.long <- rbind(pdat[,c("ID","Group","Degree")],
                   pdat[,c("ID","Group","Degree")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["irt.em.norm"]] <- c(pdat[["pre.irt.em.norm"]], pdat[["pos.irt.em.norm"]])

y.position.min <- abs(
  max(pdat.long[["irt.em.norm"]])
  - min(pdat.long[["irt.em.norm"]]))/20

lvars = as.list(c("dif.irt.em.norm","pos.irt.em.norm","pre.irt.em.norm"))
names(lvars) = unlist(lvars)

Pre-test and Post-test PairWise comparisons for: irt.em.norm ~ Group:Degree

pwc.long <- group_by(pdat.long, Group:Degree) %>%
  pairwise_wilcox_test(irt.em.norm ~ time, detailed = T)

df <- pwc.long[,c(".y.","Group:Degree","group1","group2","n1","n2","estimate",
                  "statistic","p.adj","p.adj.signif")]
.y. Group:Degree group1 group2 n1 n2 estimate statistic p.adj p.adj.signif
irt.em.norm Ctr.A:4a pre pos 5 5 1.2639448 15.0 0.675 ns
irt.em.norm Ctr.A:5a pre pos 7 7 11.4515997 40.0 0.053 ns
irt.em.norm Ctr.B:3a pre pos 28 28 0.9883975 402.0 0.876 ns
irt.em.norm Ctr.B:4a pre pos 21 21 -24.7571802 98.0 0.002 **
irt.em.norm Ctr.B:5a pre pos 19 19 -24.7571682 94.0 0.011 *
irt.em.norm Exp:3a pre pos 19 19 -10.4049400 159.0 0.539 ns
irt.em.norm Exp:4a pre pos 15 15 -28.9344905 52.5 0.012 *
irt.em.norm Exp:5a pre pos 8 8 -24.7571766 21.0 0.265 ns 
stat.test <- pwc.long %>% add_xy_position(x = "Group:Degree", fun = "max")
sidx = which(stat.test$p.adj.signif != "ns")
stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))

pdat.long[[paste0(c("Group","Degree"), collapse = ":")]] = apply(
  pdat.long[, c("Group","Degree")], 1, paste0, collapse = ":")

ggboxplot(pdat.long, x = "Group:Degree", y = "irt.em.norm",
          palette = color$prepost, fill = "time") +
  stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = T,
                     label = "{ p.adj } ({ p.adj.signif })") + xlab("")

pwc.long <- group_by(pdat.long, Group, Degree) %>%
  pairwise_wilcox_test(irt.em.norm ~ time, detailed = T)

stat.test <- pwc.long %>% add_xy_position(x = "time", fun = "mean_se")
sidx = which(stat.test$p.adj.signif != "ns")
stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))

gg <- ggline(
  pdat.long, x = "time", y = "irt.em.norm",
  color = "Degree", linetype = "Degree", shape = "Degree", size = 1.5,
  facet.by = "Group", add = c("mean_se"),
  position = position_dodge(width = 0.3), palette = color[["Degree"]])

pdat.long$xj = jitter(as.numeric(pdat.long[["time"]]), amount=.1)
pdat.long$yj = jitter(pdat.long[["irt.em.norm"]], amount = .01)

gg + geom_point(
  data = pdat.long, aes_string(x="xj",y="yj",colour="Degree"), size=0.5) +
  stat_pvalue_manual(
    stat.test, tip.length = 0, hide.ns = T, label.size = 5,
    position = position_dodge(width = 0.3), color = "Degree",
    label = "{ p.adj } ({ p.adj.signif })") + xlab("") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

Scheirer and Wilcoxon PairWise comparisons for: irt.em.norm ~ Group:Degree

sch <- lapply(lvars, FUN = function(x) {
  scheirer.test(pdat, x, c("Group","Degree"), as.table = T) 
})
df <- do.call(rbind.fill, sch)
var Effect Df Sum Sq H p.value p.value.signif
dif.irt.em.norm Group 2 17343.0632 13.9278137 0.0009454 ***
dif.irt.em.norm Degree 2 21003.9892 16.8678188 0.0002174 ***
dif.irt.em.norm Group:Degree 3 925.6556 0.7433726 0.8629579 ns
dif.irt.em.norm Residuals 114 118072.4111
pos.irt.em.norm Group 2 10204.9967 8.6280586 0.0133795 *
pos.irt.em.norm Degree 2 25121.6060 21.2396627 0.0000244 ****
pos.irt.em.norm Group:Degree 3 1102.7887 0.9323791 0.8176078 ns
pos.irt.em.norm Residuals 114 112764.7928
pre.irt.em.norm Group 2 499.8375 0.4160614 0.8121821 ns
pre.irt.em.norm Degree 2 3596.6117 2.9937950 0.2238235 ns
pre.irt.em.norm Group:Degree 3 1385.9006 1.1536142 0.7641497 ns
pre.irt.em.norm Residuals 114 138667.8594 
pwc <- lapply(lvars, FUN = function(x) {
  list(
    Group = tryCatch(pairwise_wilcox_test(group_by(pdat, Degree),
                                 as.formula(paste0(x," ~ Group")))
                         , error = function(e) NULL),
    Degree = tryCatch(pairwise_wilcox_test(group_by(pdat, Group),
                                 as.formula(paste0(x," ~ Degree")))
                         , error = function(e) NULL)
  )
})

df <- do.call(rbind.fill, lapply(pwc, FUN =  function(x) {
  do.call(rbind.fill, x)
}))

ivs = c()
if ("Group" %in% colnames(df)) ivs = c(ivs, "Group")
if ("Degree" %in% colnames(df)) ivs = c(ivs, "Degree")
df <- df[,c(".y.",ivs,"group1","group2","n1","n2",
            "statistic","p.adj","p.adj.signif")]
.y. Group Degree group1 group2 n1 n2 statistic p.adj p.adj.signif
dif.irt.em.norm 3a Ctr.B Exp 28 19 224.0 0.368 ns
dif.irt.em.norm 4a Ctr.A Ctr.B 5 21 25.5 0.240 ns
dif.irt.em.norm 4a Ctr.A Exp 5 15 20.0 0.274 ns
dif.irt.em.norm 4a Ctr.B Exp 21 15 166.5 0.783 ns
dif.irt.em.norm 5a Ctr.A Ctr.B 7 19 21.5 0.030 *
dif.irt.em.norm 5a Ctr.A Exp 7 8 8.0 0.048 *
dif.irt.em.norm 5a Ctr.B Exp 19 8 79.5 0.873 ns
dif.irt.em.norm Ctr.A 4a 5a 5 7 22.0 0.530 ns
dif.irt.em.norm Ctr.B 3a 4a 28 21 134.0 0.004 **
dif.irt.em.norm Ctr.B 3a 5a 28 19 161.5 0.048 *
dif.irt.em.norm Ctr.B 4a 5a 21 19 216.5 0.653 ns
dif.irt.em.norm Exp 3a 4a 19 15 83.0 0.121 ns
dif.irt.em.norm Exp 3a 5a 19 8 47.0 0.258 ns
dif.irt.em.norm Exp 4a 5a 15 8 64.5 0.796 ns
pos.irt.em.norm 3a Ctr.B Exp 28 19 224.5 0.369 ns
pos.irt.em.norm 4a Ctr.A Ctr.B 5 21 22.0 0.109 ns
pos.irt.em.norm 4a Ctr.A Exp 5 15 17.5 0.137 ns
pos.irt.em.norm 4a Ctr.B Exp 21 15 153.5 0.900 ns
pos.irt.em.norm 5a Ctr.A Ctr.B 7 19 26.5 0.052 ns
pos.irt.em.norm 5a Ctr.A Exp 7 8 15.5 0.322 ns
pos.irt.em.norm 5a Ctr.B Exp 19 8 86.5 0.559 ns
pos.irt.em.norm Ctr.A 4a 5a 5 7 16.0 0.871 ns
pos.irt.em.norm Ctr.B 3a 4a 28 21 126.5 0.002 **
pos.irt.em.norm Ctr.B 3a 5a 28 19 124.5 0.004 **
pos.irt.em.norm Ctr.B 4a 5a 21 19 199.0 1.000 ns
pos.irt.em.norm Exp 3a 4a 19 15 80.0 0.080 ns
pos.irt.em.norm Exp 3a 5a 19 8 54.5 0.516 ns
pos.irt.em.norm Exp 4a 5a 15 8 68.5 0.568 ns
pre.irt.em.norm 3a Ctr.B Exp 28 19 253.5 0.790 ns
pre.irt.em.norm 4a Ctr.A Ctr.B 5 21 47.5 1.000 ns
pre.irt.em.norm 4a Ctr.A Exp 5 15 35.0 1.000 ns
pre.irt.em.norm 4a Ctr.B Exp 21 15 159.5 1.000 ns
pre.irt.em.norm 5a Ctr.A Ctr.B 7 19 86.0 0.783 ns
pre.irt.em.norm 5a Ctr.A Exp 7 8 36.5 0.783 ns
pre.irt.em.norm 5a Ctr.B Exp 19 8 82.5 0.783 ns
pre.irt.em.norm Ctr.A 4a 5a 5 7 12.0 0.384 ns
pre.irt.em.norm Ctr.B 3a 4a 28 21 224.5 0.420 ns
pre.irt.em.norm Ctr.B 3a 5a 28 19 198.5 0.420 ns
pre.irt.em.norm Ctr.B 4a 5a 21 19 186.0 0.719 ns
pre.irt.em.norm Exp 3a 4a 19 15 119.0 1.000 ns
pre.irt.em.norm Exp 3a 5a 19 8 72.5 1.000 ns
pre.irt.em.norm Exp 4a 5a 15 8 61.0 1.000 ns 
plots <- lapply(lvars, FUN = function(y) {
  livs = list("Group", "Degree")
  names(livs) = unlist(livs)
  lapply(livs, FUN = function(x) {
    iv2 = setdiff(names(livs), x)
    if (!is.null(pwc[[y]][[iv2]])) {
      stat.test <- pwc[[y]][[iv2]] %>% add_xy_position(x = x, fun = "max")
      sidx = which(stat.test$p.adj.signif != "ns")
      stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))
      
      ggboxplot(pdat, x = x, y = y, fill = iv2, palette = color[[iv2]]) +
        stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = T, label.size = 5,
                           label="{ p.adj } ({ p.adj.signif })") + xlab("")
    }
  })
})
if (!is.null(plots[["pre.irt.em.norm"]][["Group"]]) &&
    !is.null(plots[["pos.irt.em.norm"]][["Group"]])) {
  egg::ggarrange(plots[["pre.irt.em.norm"]][["Group"]],
                 plots[["pos.irt.em.norm"]][["Group"]], nrow = 1)  
}

if (!is.null(plots[["pre.irt.em.norm"]][["Degree"]]) &&
    !is.null(plots[["pos.irt.em.norm"]][["Degree"]])) {
  egg::ggarrange(plots[["pre.irt.em.norm"]][["Degree"]],
                 plots[["pos.irt.em.norm"]][["Degree"]], nrow = 1)
}

psch = sch[["dif.irt.em.norm"]]
idx = which(psch$Effect == "Group:Degree") 

dof = floor(as.double(psch$Df[idx]))
dof.res = floor(as.double(psch$Df[which(psch$Effect == "Residuals")]))
statistic = round(as.double(psch$H[idx]), 3)
p = round(as.double(psch[["p.value"]][idx]), 3)
pval = ifelse(p < 0.001,paste0(" , p<0.001"),paste0(" , p=",p))

if (!is.null(plots[["dif.irt.em.norm"]][["Group"]]))
  plots[["dif.irt.em.norm"]][["Group"]] +
    labs(subtitle = paste0("Scheirer-Ray-Hare H(", dof, ",", 
          dof.res, ")=", statistic, pval),
         caption = get_pwc_label(pwc[["dif.irt.em.norm"]][["Degree"]])) +
    ylab("irt.em.norm (dif)") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

psch = sch[["dif.irt.em.norm"]]
idx = which(psch$Effect == "Group:Degree") 

dof = floor(as.double(psch$Df[idx]))
dof.res = floor(as.double(psch$Df[which(psch$Effect == "Residuals")]))
statistic = round(as.double(psch$H[idx]), 3)
p = round(as.double(psch[["p.value"]][idx]), 3)
pval = ifelse(p < 0.001,paste0(" , p<0.001"),paste0(" , p=",p))

if (!is.null(plots[["dif.irt.em.norm"]][["Degree"]]))
  plots[["dif.irt.em.norm"]][["Degree"]] +
    labs(subtitle = paste0("Scheirer-Ray-Hare H(", dof, ",", 
          dof.res, ")=", statistic, pval),
         caption = get_pwc_label(pwc[["dif.irt.em.norm"]][["Group"]])) +
    ylab("irt.em.norm (dif)") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

Two-way factor analysis for: irt.em.norm ~ Group:qtl.irt.em.norm

pdat = remove_group_data(
  dat[!is.na(dat[["Group"]]) & !is.na(dat[["qtl.irt.em.norm"]]),],
  "dif.irt.em.norm", c("Group","qtl.irt.em.norm"))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning in `stats::qt()`:
## ! NaNs produced
pdat.long <- rbind(pdat[,c("ID","Group","qtl.irt.em.norm")],
                   pdat[,c("ID","Group","qtl.irt.em.norm")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["irt.em.norm"]] <- c(pdat[["pre.irt.em.norm"]], pdat[["pos.irt.em.norm"]])

y.position.min <- abs(
  max(pdat.long[["irt.em.norm"]])
  - min(pdat.long[["irt.em.norm"]]))/20

lvars = as.list(c("dif.irt.em.norm","pos.irt.em.norm","pre.irt.em.norm"))
names(lvars) = unlist(lvars)

Pre-test and Post-test PairWise comparisons for: irt.em.norm ~ Group:qtl.irt.em.norm

pwc.long <- group_by(pdat.long, Group:qtl.irt.em.norm) %>%
  pairwise_wilcox_test(irt.em.norm ~ time, detailed = T)

df <- pwc.long[,c(".y.","Group:qtl.irt.em.norm","group1","group2","n1","n2","estimate",
                  "statistic","p.adj","p.adj.signif")]
.y. Group:qtl.irt.em.norm group1 group2 n1 n2 estimate statistic p.adj p.adj.signif
irt.em.norm Ctr.A:2nd pre pos 5 5 22.365413 23.0 0.036 *
irt.em.norm Ctr.A:3rd pre pos 8 8 4.384584 40.0 0.399 ns
irt.em.norm Ctr.B:1st pre pos 12 12 -10.404944 60.0 0.498 ns
irt.em.norm Ctr.B:2nd pre pos 33 33 -29.199056 296.0 0.001 **
irt.em.norm Ctr.B:3rd pre pos 23 23 4.441855 276.0 0.796 ns
irt.em.norm Exp:1st pre pos 9 9 -10.404991 31.0 0.418 ns
irt.em.norm Exp:2nd pre pos 21 21 -26.003887 131.5 0.025 *
irt.em.norm Exp:3rd pre pos 12 12 -24.757073 49.0 0.166 ns 
stat.test <- pwc.long %>% add_xy_position(x = "Group:qtl.irt.em.norm", fun = "max")
sidx = which(stat.test$p.adj.signif != "ns")
stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))

pdat.long[[paste0(c("Group","qtl.irt.em.norm"), collapse = ":")]] = apply(
  pdat.long[, c("Group","qtl.irt.em.norm")], 1, paste0, collapse = ":")

ggboxplot(pdat.long, x = "Group:qtl.irt.em.norm", y = "irt.em.norm",
          palette = color$prepost, fill = "time") +
  stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = T,
                     label = "{ p.adj } ({ p.adj.signif })") + xlab("")

pwc.long <- group_by(pdat.long, Group, qtl.irt.em.norm) %>%
  pairwise_wilcox_test(irt.em.norm ~ time, detailed = T)

stat.test <- pwc.long %>% add_xy_position(x = "time", fun = "mean_se")
sidx = which(stat.test$p.adj.signif != "ns")
stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))

gg <- ggline(
  pdat.long, x = "time", y = "irt.em.norm",
  color = "qtl.irt.em.norm", linetype = "qtl.irt.em.norm", shape = "qtl.irt.em.norm", size = 1.5,
  facet.by = "Group", add = c("mean_se"),
  position = position_dodge(width = 0.3), palette = color[["qtl.irt.em.norm"]])

pdat.long$xj = jitter(as.numeric(pdat.long[["time"]]), amount=.1)
pdat.long$yj = jitter(pdat.long[["irt.em.norm"]], amount = .01)

gg + geom_point(
  data = pdat.long, aes_string(x="xj",y="yj",colour="qtl.irt.em.norm"), size=0.5) +
  stat_pvalue_manual(
    stat.test, tip.length = 0, hide.ns = T, label.size = 5,
    position = position_dodge(width = 0.3), color = "qtl.irt.em.norm",
    label = "{ p.adj } ({ p.adj.signif })") + xlab("") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

Scheirer and Wilcoxon PairWise comparisons for: irt.em.norm ~ Group:qtl.irt.em.norm

sch <- lapply(lvars, FUN = function(x) {
  scheirer.test(pdat, x, c("Group","qtl.irt.em.norm"), as.table = T) 
})
df <- do.call(rbind.fill, sch)
var Effect Df Sum Sq H p.value p.value.signif
dif.irt.em.norm Group 2 8017.7590 6.3425879 0.0419493 *
dif.irt.em.norm qtl.irt.em.norm 2 6574.1092 5.2005636 0.0742526 ns
dif.irt.em.norm Group:qtl.irt.em.norm 3 6703.0291 5.3025480 0.1509371 ns
dif.irt.em.norm Residuals 115 130115.3448
pos.irt.em.norm Group 2 5396.9226 4.5151228 0.1046053 ns
pos.irt.em.norm qtl.irt.em.norm 2 27899.2532 23.3408118 0.0000085 ****
pos.irt.em.norm Group:qtl.irt.em.norm 3 4269.4145 3.5718376 0.3115639 ns
pos.irt.em.norm Residuals 115 112149.3954
pre.irt.em.norm Group 2 219.7309 0.1807314 0.9135970 ns
pre.irt.em.norm qtl.irt.em.norm 2 125574.1638 103.2863287 0.0000000 ****
pre.irt.em.norm Group:qtl.irt.em.norm 3 345.6199 0.2842767 0.9629580 ns
pre.irt.em.norm Residuals 115 17200.6493 
pwc <- lapply(lvars, FUN = function(x) {
  list(
    Group = tryCatch(pairwise_wilcox_test(group_by(pdat, qtl.irt.em.norm),
                                 as.formula(paste0(x," ~ Group")))
                         , error = function(e) NULL),
    qtl.irt.em.norm = tryCatch(pairwise_wilcox_test(group_by(pdat, Group),
                                 as.formula(paste0(x," ~ qtl.irt.em.norm")))
                         , error = function(e) NULL)
  )
})

df <- do.call(rbind.fill, lapply(pwc, FUN =  function(x) {
  do.call(rbind.fill, x)
}))

ivs = c()
if ("Group" %in% colnames(df)) ivs = c(ivs, "Group")
if ("qtl.irt.em.norm" %in% colnames(df)) ivs = c(ivs, "qtl.irt.em.norm")
df <- df[,c(".y.",ivs,"group1","group2","n1","n2",
            "statistic","p.adj","p.adj.signif")]
.y. Group qtl.irt.em.norm group1 group2 n1 n2 statistic p.adj p.adj.signif
dif.irt.em.norm 1st Ctr.B Exp 12 9 54.0 1.00e+00 ns
dif.irt.em.norm 2nd Ctr.A Ctr.B 5 33 18.0 1.50e-02 *
dif.irt.em.norm 2nd Ctr.A Exp 5 21 9.0 1.50e-02 *
dif.irt.em.norm 2nd Ctr.B Exp 33 21 383.5 5.17e-01 ns
dif.irt.em.norm 3rd Ctr.A Ctr.B 8 23 85.0 1.00e+00 ns
dif.irt.em.norm 3rd Ctr.A Exp 8 12 38.5 1.00e+00 ns
dif.irt.em.norm 3rd Ctr.B Exp 23 12 124.5 1.00e+00 ns
dif.irt.em.norm Ctr.A 2nd 3rd 5 8 9.0 1.22e-01 ns
dif.irt.em.norm Ctr.B 1st 2nd 12 33 177.0 5.98e-01 ns
dif.irt.em.norm Ctr.B 1st 3rd 12 23 170.0 5.30e-01 ns
dif.irt.em.norm Ctr.B 2nd 3rd 33 23 563.0 7.00e-03 **
dif.irt.em.norm Exp 1st 2nd 9 21 88.0 1.00e+00 ns
dif.irt.em.norm Exp 1st 3rd 9 12 57.0 1.00e+00 ns
dif.irt.em.norm Exp 2nd 3rd 21 12 157.0 7.53e-01 ns
pos.irt.em.norm 1st Ctr.B Exp 12 9 57.0 8.54e-01 ns
pos.irt.em.norm 2nd Ctr.A Ctr.B 5 33 22.0 1.50e-02 *
pos.irt.em.norm 2nd Ctr.A Exp 5 21 9.0 1.30e-02 *
pos.irt.em.norm 2nd Ctr.B Exp 33 21 332.0 7.96e-01 ns
pos.irt.em.norm 3rd Ctr.A Ctr.B 8 23 83.5 1.00e+00 ns
pos.irt.em.norm 3rd Ctr.A Exp 8 12 38.5 1.00e+00 ns
pos.irt.em.norm 3rd Ctr.B Exp 23 12 126.5 1.00e+00 ns
pos.irt.em.norm Ctr.A 2nd 3rd 5 8 0.5 5.00e-03 **
pos.irt.em.norm Ctr.B 1st 2nd 12 33 98.5 1.80e-02 *
pos.irt.em.norm Ctr.B 1st 3rd 12 23 58.0 1.30e-02 *
pos.irt.em.norm Ctr.B 2nd 3rd 33 23 313.5 2.52e-01 ns
pos.irt.em.norm Exp 1st 2nd 9 21 41.0 2.90e-02 *
pos.irt.em.norm Exp 1st 3rd 9 12 16.0 1.80e-02 *
pos.irt.em.norm Exp 2nd 3rd 21 12 98.0 2.74e-01 ns
pre.irt.em.norm 1st Ctr.B Exp 12 9 69.5 2.31e-01 ns
pre.irt.em.norm 2nd Ctr.A Ctr.B 5 33 81.0 9.74e-01 ns
pre.irt.em.norm 2nd Ctr.A Exp 5 21 41.5 9.74e-01 ns
pre.irt.em.norm 2nd Ctr.B Exp 33 21 275.0 5.91e-01 ns
pre.irt.em.norm 3rd Ctr.A Ctr.B 8 23 88.0 1.00e+00 ns
pre.irt.em.norm 3rd Ctr.A Exp 8 12 44.0 1.00e+00 ns
pre.irt.em.norm 3rd Ctr.B Exp 23 12 133.0 1.00e+00 ns
pre.irt.em.norm Ctr.A 2nd 3rd 5 8 0.0 1.00e-03 **
pre.irt.em.norm Ctr.B 1st 2nd 12 33 0.0 3.00e-07 ****
pre.irt.em.norm Ctr.B 1st 3rd 12 23 0.0 1.00e-07 ****
pre.irt.em.norm Ctr.B 2nd 3rd 33 23 0.0 0.00e+00 ****
pre.irt.em.norm Exp 1st 2nd 9 21 0.0 3.64e-05 ****
pre.irt.em.norm Exp 1st 3rd 9 12 0.0 3.64e-05 ****
pre.irt.em.norm Exp 2nd 3rd 21 12 0.0 4.80e-06 **** 
plots <- lapply(lvars, FUN = function(y) {
  livs = list("Group", "qtl.irt.em.norm")
  names(livs) = unlist(livs)
  lapply(livs, FUN = function(x) {
    iv2 = setdiff(names(livs), x)
    if (!is.null(pwc[[y]][[iv2]])) {
      stat.test <- pwc[[y]][[iv2]] %>% add_xy_position(x = x, fun = "max")
      sidx = which(stat.test$p.adj.signif != "ns")
      stat.test$y.position[sidx] <- stat.test$y.position[sidx] + y.position.min * (1:length(sidx))
      
      ggboxplot(pdat, x = x, y = y, fill = iv2, palette = color[[iv2]]) +
        stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = T, label.size = 5,
                           label="{ p.adj } ({ p.adj.signif })") + xlab("")
    }
  })
})
if (!is.null(plots[["pre.irt.em.norm"]][["Group"]]) &&
    !is.null(plots[["pos.irt.em.norm"]][["Group"]])) {
  egg::ggarrange(plots[["pre.irt.em.norm"]][["Group"]],
                 plots[["pos.irt.em.norm"]][["Group"]], nrow = 1)  
}

if (!is.null(plots[["pre.irt.em.norm"]][["qtl.irt.em.norm"]]) &&
    !is.null(plots[["pos.irt.em.norm"]][["qtl.irt.em.norm"]])) {
  egg::ggarrange(plots[["pre.irt.em.norm"]][["qtl.irt.em.norm"]],
                 plots[["pos.irt.em.norm"]][["qtl.irt.em.norm"]], nrow = 1)
}

psch = sch[["dif.irt.em.norm"]]
idx = which(psch$Effect == "Group:qtl.irt.em.norm") 

dof = floor(as.double(psch$Df[idx]))
dof.res = floor(as.double(psch$Df[which(psch$Effect == "Residuals")]))
statistic = round(as.double(psch$H[idx]), 3)
p = round(as.double(psch[["p.value"]][idx]), 3)
pval = ifelse(p < 0.001,paste0(" , p<0.001"),paste0(" , p=",p))

if (!is.null(plots[["dif.irt.em.norm"]][["Group"]]))
  plots[["dif.irt.em.norm"]][["Group"]] +
    labs(subtitle = paste0("Scheirer-Ray-Hare H(", dof, ",", 
          dof.res, ")=", statistic, pval),
         caption = get_pwc_label(pwc[["dif.irt.em.norm"]][["qtl.irt.em.norm"]])) +
    ylab("irt.em.norm (dif)") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))

psch = sch[["dif.irt.em.norm"]]
idx = which(psch$Effect == "Group:qtl.irt.em.norm") 

dof = floor(as.double(psch$Df[idx]))
dof.res = floor(as.double(psch$Df[which(psch$Effect == "Residuals")]))
statistic = round(as.double(psch$H[idx]), 3)
p = round(as.double(psch[["p.value"]][idx]), 3)
pval = ifelse(p < 0.001,paste0(" , p<0.001"),paste0(" , p=",p))

if (!is.null(plots[["dif.irt.em.norm"]][["qtl.irt.em.norm"]]))
  plots[["dif.irt.em.norm"]][["qtl.irt.em.norm"]] +
    labs(subtitle = paste0("Scheirer-Ray-Hare H(", dof, ",", 
          dof.res, ")=", statistic, pval),
         caption = get_pwc_label(pwc[["dif.irt.em.norm"]][["Group"]])) +
    ylab("irt.em.norm (dif)") +
  theme(strip.text = element_text(size = 16),
        axis.text = element_text(size = 18))